#### Regressions that compare effects on yield, produce, and harvest

## Setup environment

setwd("~/research/coffee-dataverse/src")

do.origspec <- T
source("load.R", encoding="iso-8859-1")

library(lfe)
library(splines)

## Prepare data

df$logharvest <- log(df$total.harvest)
df$logharvest[!is.finite(df$logharvest)] <- NA

df$logproduce <- log(df$total.produce)
df$logproduce[!is.finite(df$logproduce)] <- NA

## Run main regressions

formula.rhs <- paste("~", mainspec.rhs, "| 0 | Munic�pio + year")
mod.yield <- felm(as.formula(paste("logyield", formula.rhs)), data=df, keepCX=T)
mod.produce <- felm(as.formula(paste("logproduce", formula.rhs)), data=df, keepCX=T)
mod.harvest <- felm(as.formula(paste("logharvest", formula.rhs)), data=df, keepCX=T)

## Print the regression tables

library(stargazer)

stargazer(list(mod.yield, mod.produce, mod.harvest), add.lines=list(c('Municipality FE', rep('Yes', 3)),
                                                       c('Year FE', rep('Yes', 3)),
                                                       c('State Trends', rep('Cubic', 3))),
          dep.var.labels=c("Log Yields", "Log Production", "Log Harvests"), covariate.labels=weatherlabels,
          omit=':', df=F) #, out="../tables/prodharvmods.tex")

source("conley.R")

mod.yield.cse <- ConleySEs(mod.yield, df[-mod.yield$na.action,], "Munic�pio", "year", "latitude", "longitude", dist_cutoff=400, lag_cutoff=1)
mod.produce.cse <- ConleySEs(mod.produce, df[-mod.produce$na.action,], "Munic�pio", "year", "latitude", "longitude", dist_cutoff=400, lag_cutoff=1)
mod.harvest.cse <- ConleySEs(mod.harvest, df[-mod.harvest$na.action,], "Munic�pio", "year", "latitude", "longitude", dist_cutoff=400, lag_cutoff=1)

stargazer(list(mod.yield, mod.produce, mod.harvest), add.lines=list(c('Municipality FE', rep('Yes', 3)),
                                                       c('Year FE', rep('Yes', 3)),
                                                       c('State Trends', rep('Cubic', 3))),
          dep.var.labels=c("Log Yields", "Log Production", "Log Harvests"), covariate.labels=weatherlabels,
          omit=':', df=F,
          se=list(sqrt(diag(mod.yield.cse$Spatial_HAC)), sqrt(diag(mod.produce.cse$Spatial_HAC)), sqrt(diag(mod.harvest.cse$Spatial_HAC))))
